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Einstein condensates. Among these families, two are always stable, while the third one is only stable 
sufficiently close to the linear regime. The solitons' excitation spectra reveal the potential existence 
of a second anomalous mode. While the first such mode describes the soliton oscillatory motion 
in a parabolic trap, the second, when present, reflects the double well structure of the underlying 
single-particle spectrum. This novel mode results in moving density stripes in the vicinity of the 
soliton core, or in an out-of-phase oscillation of the constituent components, with little effect on the 
nearly stationary striped total density of the composite soliton. 



INTRODUCTION 

The recent experimental realization of spin-orbit coupling (SOC) in Bose-Einstein condensates (BECs) [l|-|3[ and 
fermionic gases [J, l5| has stimulated considerable research interest. This is due to the fact that it paves the way 
towards a deeper understanding of exotic solid state systems such as topological insulators, making use of the highly 
controllable environment of cold atom physics [a, lZ| . Ultracold atomic gases promise to be an ideal system for the 
generation of versatile artificial gauge fields [8[ and studying such Abelian or even non-Abclian gauge fields is relevant 
not only in mimicking solid state physics, but potentially even to fundamental theories such as quantum electro- or 
chromodynamics [9[. 

In the above context, ground state properties of SOC-BECs, inclu ding phase separation and the existence of the 



stripe phase [10|, as well as the collective dynamics and oscillations [llj, have already been studied in some detail 
both in theory and experiments [l|, y]- Furthermore, a relevant direction that has also attracted much attention is the 
study of topological excitations, such as skyrmions [12|, vortices [13|, and Dirac monopoles IJ]; additionally, dark 
solitons in toroidal geometry [15[ and bright solitons in quasi one-dimensional (ID) attractive SOC-BECs [16[ have 
also been predicted. While the above structures have already been studied extensively in the context of single- and 



multi-component BECs |17lll8l|. the presence of SOC can significantly enrich their structural, stability and dynamical 
properties. 

It is the purpose of this work to highlight the above by presenting and analyzing dark soliton states in SOC-BECs 
confined in a highly anisotropic (quasi-lD) parabolic trap. Employing a multiscale expansion method, we find three 
different dark soliton families, featuring either a constant or a spatially modulated background density; in the latter 
case, dark solitons occur as excited states on top of the stripe phase of SOC-BECs (we call these states "stripe 
solitons"). We perform a Bogolyubov-de Gennes (BdG) linearization analysis of the above soliton families, showing 
that constant background solitons are always stable, while stripe solitons are stable only close to the linear limit. The 
characteristic negative energy (anomalous) mode associated with the oscillation frequency of solitons near the trap 
center is identified within the excitation spectrum and its eigenfrequency is determined analytically. Importantly, in 
certain parameter regions we also find a second anomalous mode, which does not exist in single- or multi-component 
BECs, and is only sustained due to the double well structure of the SOC single particle spectrum, featuring two 
minima. Exciting this mode, we observe intriguing dynamics: constant background solitons feature moving stripes in 
the vicinity of the soliton core, somewhat reminiscent of the Kelvin mode modulation of vortex lines [19[; on the other 
hand, in the case of stripe solitons, we observe an out-of-phase oscillation of the constituent solitons, which does not 
affect the stationary total density of the composite state. 

THE MODEL AND ITS ANALYTICAL CONSIDERATION 

We consider a SOC-BEC confined in a quasi-lD parabolic trap, with longitudinal and transverse frequencies u)x <C 
uj± . In the framework of mean- field theory, this system can be described by the energy functional [l|, [lOj : 

E^n^UoM + \ {gxx\uf -f g^M^ + 2grMM^) , (1) 

where u = {u, v)'^ , and the condensate wavefunctions u and v are related (through suitable rotations [l|) to the two 
pseudo-spin components of the BEC. The single particle Hamiltonian Hq in Eq. ([ij reads: Hq = ^{Px^ + kLO'z)'^ + 
Vti(x)l -|- flax + ScTz, where Px = —ihdx is the momentum operator in the longitudinal direction, m is the atomic 
mass, ax,z are the Pauli matrices, 1 is the unit matrix, k^ is the wavenumber of the Raman laser which couples the 
states u and v, S is the detuning from Raman resonance, and fi is the strength of the Raman coupling. Additionally, 
Vtr(a::) is the external trapping potential, considered to assume the usual parabolic form: Vtr = (l/2)r7T,a;^a;^. Finally, 
the effectively ID coupling constants g^, arc given by gij = 2hLL)±aij, where aij are the s-wave scattering lengths 
(assumed to be positive). Measuring energy in units of huj±, length in units of the transverse harmonic oscillator 
length a± = y^h/{mLL!±), time in units of wj , and densities in units of the scattering length an, we derive from 
Eq. ^ the following dimensionless equations of motion: 

idtu = (~^dl - ikLdx + Vt, + |up + ,5|if + &\u + nv, 
1. 

(2) 



idtv = ( -l-dl + ikLdx + Vu + P\u\^ + l\v\^ 'S\v + nu, 



where (3 = ai2/Q;ii, 7 ~ 0^22/ olii, and wc have used k^ — >■ k^/a^, VL — > f7?ia;^ and (5 — >■ (5fia;j^; the trapping potential in 
Eqs. ([2]) is now given by \4r(a;) = (l/2)a;tra;^, where Wtr = Wa;/w_L. The stationary counterpart of Eqs. (I2|) is obtained 
by factorizing u(a;,i) = u(a;) exp(— i//i), where ^ denotes the chemical potential. 

In the following, we adopt experimentally relevant values for the above parameters. Considering, as in the case of 
the experiments of Ref. [l[, two spin states of the ^^Rb 5Si/2, F = 1 manifold, we will use the following dimensionless 
parameter values: Raman wavenumber k^ = 8, normalized Raman coupling strength fi G [0, 100], and ratios of the 
scattering lengths an : a\2 '■ C122 = 1 '■ 0.995 : 0.995. It is thus physically relevant to use below the approximation 
7^1 while we will let /3 be a free parameter. 

We now use a multiscale perturbation method to derive approximate dark soliton solutions of Eqs. ([2]). First we 
introduce the ansatz: u = [[/, y] exp[i(fcx — /ii)] , where k is the momentum and /i = w + e^ojo the chemical potential, 
where w is the energy in the linear limit, e^wo is a small deviation about this energy (0 < e ^ 1), and wo/a; = 0{\). 
We also assume that the trap strength is Wtr = e^(Dtr- Next, introducing the slow variables T — e^i, X = ex, and 
expanding f/ and ^ as f/ = eC7i(X,T) + e'^UziX.T) + . . . and ^ = eVi{X,T) + e'^V2iX,T) + . . ., we derive from 
Eqs. ([2]) the following equations at the orders 0{e^), C(e^) and 0{e^), respectively: 

Wui = 0, (3) 

Wu2 =iWoaxUi, (4) 

Wu3 = iWo^A-Ua - (idT + -d\ - A + wo j ui, (5) 

where u,; = [t/^, V^]"^ (i g {1, 2, 3}) are unknown vectors, and matrices W and A are given by: 

_ w - fcV2 - fcfci - 8 -Vl 

~Vl u-k'^l2 + kkL + i 

|C/iP + /3|Fi|2 + 14, 

/3|[/ip + |Fi|2 + -i4, 

where Wo — (W — wl)', primes denote differentiation with respect to fc, and Vtr(^) = (l/2)ii'^r^^. At ©(e^), the 
solvability condition dctW = yields the single-particle spectrum (dispersion relation) 

1. 



u = uj±{k) = -e ± ^{kkL + 5f + rj2, (6) 

displayed in Fig[T] Let L = [1, Q] and R = [1, Q]"^ be the left and right eigenvectors of W at eigenvalue 0, where 

Q^Q{u:,k) = ^u:-h^~kkL-6\. (7) 

Then, the compatibility condition of Ec^. ([3]) yields: 

ui=Ri^(X,T), (8) 

where 'ip{X,T) is an unknown scalar field. 

Next, proceeding with Eq. (U), the compatibility condition LWqR = at the order ©(e^) enforces a vanishing 
group velocity Vg = duj/dk, i.e.: 

v,^k-kL^^^O. (9) 

Thus, the compatibility condition at 0{e^) requires that the parameters u and k are evaluated at stationary points 
of the dispersion relation (|6]). In the following we focus on the energy minima (wmin, ^min)- Additionally, at the order 
0{e^), we obtain the following solution for U2: 

U2 = -iR'dx^{X,T). (10) 

Finally, at 0{e^), the compatibility condition for Eq. ([5]), combined with Eqs. ([5]) and the above result for U2, yields 
the following nonlinear Schrodinger (NLS) equation for the scalar field ip: 

idri^ + \d\^ - (y\,\,f - ujo)^P = 14r(A)^, (11) 
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FIG. 1. (Color online) Sketches of the linear energy spectrum lj = cj±(fc) for 5 = (top panels) and (5 7^ (bottom panels), 
and also for for fl > kj^ (left panels) and for $1 < fc|_ (right panels). In the bottom panels, dashed lines depict - for comparison 
- the respective curves corresponding to 5 = 0. 
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FIG. 2. (Color online) In each of the three triplets of panels, shown are the real (top left) and imaginary (bottom left) parts of 
u and V, as well as total density profiles |up + \v\^ (right) of stationary kinks for 5 — 0, jS = 1 and e^ujo = 0.2. The top triplet 
shows a soliton in Regime I for Q/ki^ = 1.4, while the middle and bottom triplets show a fcmin-soliton and a stripe soliton in 
Regime II for fi/fcf, — 0.625. Solid (black) lines depict numerical results, while dashed (green and yellow) lines correspond to 
analytical results [cf. Eqs. p3l) and psp ]. 



where the coefficients A and v defined as: 



A = l 



2QQ'(fcL - k) 



1 + Q2 



(12) 



are evaluated at the minima (w, k) ~ (wmin, fcmin) of the dispersion relation ^. 

In the homogeneous case (wtr = 0), for Xv > 0, and assuming the boundary conditions \'ip\'^ — > LOQ/iy as \X\ — >■ 00, 
the above scalar NLS Eq. (|TT]) possesses dark soliton solutions, ip-os, characterized by the free parameter luq. A single 
dark soliton reads: 



"^DS = \/ ^qJv [cos 0tanh(?7) + i sin S\ , 



(13) 



where j] = \Ju)q/\ cos 9[X — Xq{T)], and 9 is the so-called "soliton phase angle" controlling the velocity and the phase 
shift across the dark soliton {\9\ < tt/2). Additionally, Xo{T) is the soliton center, while the soliton amplitude (depth) 
and soliton velocity are respectively given by -^/wo/Acos^ and Xq — \/u}oJXsm9. Note that the limiting case 6* = 
corresponds to a stationary kink ("black" soliton), while cases with 0^0 give rise to traveling ("grey") solitons. 
Using the above expressions, we can now write down an approximate dark soliton solution of Eqs. ^, in terms of 



the original variables x and t, as follows: 

f j w e-0Dsexp [ifcmina;-i(a;min + e^cjo)t]fg j, (14) 

where Tp^c, is given in Eq. p^ . but with argument 77 -^ rj/e, and Q„un = Q{uJmim ^min)- 

DARK SOLITONS IN THE HOMOGENEOUS SETTING 

The coefficients of the NLS equation (jlip and thus the soliton parameters, explicitly depend on Qmin, which is 
calculated at the minimum of the dispersion relation ([6]) . The latter possesses an upper branch w+ and a lower branch 
id-, as shown in Fig. [T] We hereafter focus our analysis around the energy minima of the lower branch w_. We will 
study separately the regimes fl/fc|, > 1 (Regime I) and r2/fc| < 1 (Regime II), for both 6 = and 6^0, which feature 
different characteristics regarding the energy minima, as shown in the left and right panels of Fig. [l] respectively. 

Dark solitons in Regime I. In this regime, and ioi 6 = 0, the lower branch possesses a minimum uj^^J^-^ = — J7 at zero 
momentum fcmin = 0, while Qmin = —1 and Q'^^^ = — /cl/^- The above values determine A and v in Eqs. ([T^ and, 
thus, the form of the soliton in Eq. ([T5|) . The existence of the relevant soliton solutions is numerically confirmed by 
solving the stationary version of Eqs. ([2]), where fi = — fl + e^wo [as per Eq. ((T3)) ]. using a fixed-point algorithm [20|. A 
typical example of a stationary kink {6 = 0), corresponding to parameter values ri/fc| = 1.4 and e^wp = 0.2, is shown 
in the top triplet of panels of Fig. [2j It is observed that the real parts of u and v (top left panel of the triplet) are 
opposite, in accordance with the form of the right eigenvector R = [1,-1], the imaginary parts (bottom left panel 
of the triplet) are of order 0{€^), and are equal having a sech^ profile, in accordance to the analytical prediction for 
U2. Finally, the spatial profile of the total density |up + |i;p (top right panel) has the form of a scalar dark soliton's 
density ^. 

Dark solitons in Regime II. In this regime, and for (5 = 0, we find solitons with energies near the energy minimum 
'^min ^t finite momenta, namely fcmin = ±fcL(l — f7^/fc|)^/^ (cf. top right panel of Fig. [T]). For these values of the 
energy and momentum, Q^^l^ = ~k^in^~-^{kL ± fcmin) and Q,jji„ = — r2~^(fci ± fcmin), with these values determining 
the soliton parameters in Eq. (fT2|) . 

It is thus clear that, in this regime, two different soliton solutions can be found, each corresponding to the locations 
k = ifcmin of the energy minimum; these will be called hereafter zbfcmin-solitons. An example of a fcmin-soliton, for 
n/fc£ = 0.625 and c'^uq = 0.2, is shown in the middle triplet of panels of Fig.[2l As can be seen, the real and imaginary 
parts of u and v have different amplitudes in this case, due to the form of R = [1, Qminl"^' furthermore, Re(w, v) and 
Im(u, v) are now spatially oscillatory with a wavelength 2'K/k^in- On the other hand, the total density profile features 
the usual (i.e., unmodulated) density dip. 

Importantly, still another branch of dark soliton solutions can also be found in Regime II as follows. First we note 
that the respective linear problem in this regime admits solutions which are linear combinations of plane waves of 
momenta k = ztfemin- Then, employing continuation arguments, we construct approximate dark soliton solutions in 
the form of a linear combination of the above mentioned zt/cmin-solitons. Solutions satisfying the symmetry u = —v 
(bar denotes complex conjugation) can be expressed as: 

, g^W; ( ^+ cos{k^inx) + iq- sin(fc,nina;) \ ,^^. 

^\-'q+ cos(fc,nina;) + iq- sin(fc,nina;)y' 

where q± = Q^^ + Q,ni„ and C is a free parameter. The existence of these solitons was also confirmed numerically, 
and a pertinent example is shown in the bottom triplet of panels of Fig. ^ for O/fcJ = 0.625, e^wo = 0.2 and C = 0.8. 
It is observed that the soliton background features a spatially modulated density, reminiscent of the stripe phase of 
the SO-coupled BECs [lO!]; for this reason, solitons of this branch will be called "stripe solitons". 

For all types of solitons found in Regimes I and II, the analytical approximations of Eqs. (fT3|l and (fT5|) are in 
an excellent agreement with the numerically obtained solutions (see, respectively, dashed and solid lines in Fig. [2]) . 
Furthermore, apart from the case of stationary (black) solitons, we have also studied traveling (grey) solitons with 
6* 7^ 0; this was done by direct numerical integration of Eqs. dH, by means of a Runge-Kutta method, using as initial 
conditions the analytical form of moving solitons, cf. Eq. (|13p . We have confirmed that near the linear limit (i.e., for 
e^wo ^0.1) the solitons remain robust and evolve without distortion (results not shown here). 

The case of nonzero detuning parameter. Dark solitons can also be found for 6 =/= 0, upon calculating the energy 
minimum (fcmin, Wmin) of the lower branch of the energy spectrum (jS]) (cf. bottom panels of Fig. [IJ and then using 
Eqs. P^ to determine the soliton parameters for the solution ([T5|) . There are two important differences between the 
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FIG. 3. (Color online) Top panel: the lowest imaginary eigenvalues of the linearization spectrum as functions of /i, for a dark 
soliton in Regime I, with 5 = 0, f2/fc|, = 1.4, and (itr = 0.1. Circles (in red) denote the anomalous mode eigenfrequency and 
the solid (black) line depicts the prediction of Eq. (|16fl . Bottom panel: contour plot showing the evolution of the density of a 
SOC-BEC carrying a dark soliton initially placed at a;o = 1. Dashed (white) line depicts the analytical result of Eq. (|16|l . 



cases (5^0 and (5 = 0. First, since fcmin ^ for every (5 7^ (see Fig. [T]), the dark solitons in Regime I are shifted 
to finite fc. Second, since in Regime II the degeneracy of the energy minima in the lower branch is always lifted for 
(5 7^ 0, stripe solitons cannot be constructed using the linear superposition argument as in Eq. ([T5|) . In both regimes, 
we have numerically confirmed the existence of dark solitons for (5 7^ at the global minimum of the energy spectrum. 
As a result, in both regimes I and II, these structures share similar profiles to the ones presented in the middle triplet 
of panels of Fig. [21 hence they will not be shown here. 



STABILITY AND DYNAMICS OF SOLITONS IN THE TRAP 



We now focus on the experimentally relevant setting where the parabolic trap is present (wtr 7^ 0)- In this 
case, we have confirmed that all families of soliton solutions presented above persist. Furthermore, for the sta- 
tionary solutions, we have performed a linear stability analysis and identified regimes of (in)stability. Our anal- 
ysis relies on the study of the BdG excitation spectrum around a stationary soliton solution Ugoi = (wsoi, Wsoi)^ 
of Eqs. ([2]), evolving at the chemical potential /i. The spectrum is obtained as follows. We introduce the ansatz 
u — {usoi -I- £[exp(Ai)a(a;) -I- exp(Ai)b(a;)]} exp(— i/ii), where e is a formal small parameter, and {A, (a, b)} define an 
eigenvalue-eigenvector pair. Then, substituting the above ansatz into Eqs. ^^ and linearizing, we arrive at 0(e) at the 
ensuing eigenvalue problem for eigenvectors (a, b) and eigenvalues A. Notice that as the latter are generally complex, 
i.e., A = Ar -I- JAi, instability corresponds to A^ > 0. 

First, we consider the case of solitons in Regime I, for (5 = 0, and study their stability starting from the linear limit, 
/i = — J7, and entering into the nonlinear regime by increasing the parameter ojo. In the spectrum of the respective 
soliton branch, no real eigenvalues appear and, thus, this branch is dynamically stable. As shown in the top panel 
of Fig. |3] for ^lk\ = 1.4 and Wtr = 0.1, there exists one anomalous (negative energy) mode, depicted by red circles, 
whose frequency characterizes the small-amplitude oscillations of the dark soliton around the trap center. The latter 
can be obtained analytically in the framework of Eq. (fTT|) as follows. As is well known, sufficiently deep (almost black) 
dark solitons oscillate in a parabolic trap of strength wtr with a frequency Wgoi = i^tr/vS [l8|. Hence, we can infer 
that solitons of Eq. (|TT|) with 6* « evolve in the trap so that their center ATq (t) satisfies the equation of motion: 
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(16) 



The above soliton oscillation frequency Wsoi is depicted by the solid (black) line in the top panel of Fig. [H It can be 
observed that this analytical result almost coincides with the anomalous mode eigenfrequency. To further elaborate on 
this result, we have numerically integrated Eqs. ([2]) with an initial condition of a soliton with = 0, initially placed at 
xo = 1. As seen in the bottom panel of Fig. [31 the soliton indeed performs harmonic oscillations accurately described 
by Eq. (fT6)) . cf. dashed (white) line in the figure. We note that, for (5 7^ 0, solitons have a spectrum qualitatively 
similar to the one in Fig. [3l featuring one anomalous mode and no real (i.e., unstable) eigenvalues. 

Next, we study the stability of solitons in Regime II (for b = 0), in the presence of the trap. First we note that 
^min-solitons are stable, as no real eigenvalues appear in the spectrum - see the example shown in the top panel of 
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FIG. 4. (Color online) Top panel: same as the top panel of Fig. [3l but for the fcmin-solitons in Regime II with il/kj^ = 0.625. 
Bottom panel: contour plot showing the evolution of the density of the u-component, when perturbed by the eigenvectors of 
the second anomalous mode at /^ = —43.5. Other parameter values are the same as in Fig. (3] 



Fig. m for a soliton with fi//s| = 0.625. However, an important difference arises in this case, namely the emergence 
of a second anomalous mode in the linearization spectrum. Indeed, as shown in the figure [cf. circles (in red)], there 
exist two anomalous modes: a lower- lying one (with smaller eigcnfrcquencies), associated with the soliton oscillations 
as before, and an upper-lying anomalous mode (with larger eigcnfrcquencies). This second anomalous mode docs not 
exist in the case of dark solitons in single- or multi-component BECs: its emergence is particular to dark solitons 
in SOC-BECs, which feature a double well structure in their energy spectrum in Regime II. Indeed, in the linear 
limit, the system possesses two eigenstates that energetically lie below the fcmin-soliton, namely one ground state 
configuration in each of the two wells of the energy spectrum. Each of these gives rise to one anomalous mode in 
the soliton's linearization spectrum. Superimposing the /cmin-soliton at momentum fcmin > with the ground state 
at fcmin > leads to the familiar soliton oscillation, represented by the lower-lying anomalous mode. On the other 
hand, superimposing the fcmin-soliton at momentum knun > with the ground state at fcmin < 0, different dynamics 
is to be expected: The addition of counter-propagating waves of opposite k will give rise to stripe signatures in the 
total density. We have numerically solved Eqs. ^, perturbed by the eigenvector (a, b) of the second anomalous 
mode, to check this prediction. The bottom panel of Fig. H] shows the corresponding evolution of the density of the 
M-componcnt (in which the stripe- forming effect of the perturbation is more evident, since it is the smaller component 
in our example). Clearly, exciting the second anomalous mode indeed induces a periodic spatial modulation of the 
background, moving also through the soliton, with a wave number ~ 2fcmin. This modulation effect inside the dark 
soliton core is reminiscent of the way in which Kelvin modes modulate vortex lines [l9| . 



Lastly, we consider stripe solitons in Regime II [cf. Eq. ([T5|)] in the presence of the trap. The BdG spectrum for 
such a solution with il/kj^ = 0.625 is depicted in Fig. [H As expected from the previous discussion, the spectrum 
also features two anomalous modes, reflecting the existence of a doubly degenerate lower energy (ground) state. The 
anomalous modes correspond to an in-phase and an out-of-phase superposition of these ground states, respectively. 
The former leads to in-phase oscillatory dynamics of the solitons in the two components. This is captured by the 
lower-lying anomalous mode in the spectrum, as we have checked. On the other hand, the second anomalous mode 
corresponds to out-of-phase oscillations of the solitons at fcmin < and fcmin > 0. This is directly confirmed by 
propagation of an initially stationary stripe soliton, perturbed by the respective eigenvector, as shown in the contours 
of Fig. [6l The dark solitons in each component are shown to perform small-amplitude oscillations; the inset in Fig. [6] 
depicts the evolution of the center of mass for the individual components, clearly revealing that these are out-of phase 
oscillations. Additionally, and contrary to zbfcmin-solitons. Fig. [5] shows that stripe solitons have a small stability 
domain only close to the linear limit, and become unstable deeper in the nonlinear regime; the instability, though, 
is found to be induced by the background {i.e., not from the anomalous modes). Apart from subsequent collisions 
with other modes producing small complex quartets (emerging as small bubbles in the bottom panel of Fig. [5]) , the 
two anomalous modes remain purely imaginary, while different background modes generate a cascade of instabilities 
through the bifurcation of real eigenvalue pairs. 
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FIG. 5. (Color online) Top and bottom panels show, respectively, the imaginary and real parts of the eigenvalues of the 
linearization spectrum as functions of /i, for a stripe soliton in Regime II for f2/fc|, = 0.625. The stripe soliton is stable, 
i.e., Re(A) = 0, only near the linear regime {fi < —44.15) and becomes unstable deeper in the nonlinear regime. Background 
instabilities lead to real pairs, while collisions of the anomalous modes with other modes produce weak instabilities via eigenvalue 
quartets (small bubbles of the bottom panel). 




FIG. 6. (Color online) Same as in bottom panel of Fig. |4] but for both components of the stripe soliton in Regime II with 
/i = —44. The inset middle panel shows the time evolution of the center of mass Xcm of the u component (initially at Xcm > 
- cf. blue line), of the v component (initially at Xcm < - cf. red line), and of the total density (black line), at Xcm — 0. 



CONCLUSIONS 

In summary, we have studied the existence, stability and dynamics of dark solitons in spin-orbit coupled BECs. We 
developed a perturbative approach to obtain solitons on top of either constant or spatially modulated background 
density (stripe solitons). A linear stability analysis has shown that constant background solitons are always stable, 
while stripe solitons are stable only close to the linear limit. The eigenfrequency of the anomalous mode associated with 
the oscillatory motion of solitons in a parabolic trap was determined analytically. Importantly, in a certain parameter 
regime where the single particle spectrum features a double well structure, we found a second anomalous mode, which 
does not exist in single- or multi-component BECs. Exciting this mode, we found that constant background soHtons 
feature a periodic structure of moving stripes in the vicinity of the soliton core, reminiscent of the Kelvin modulation 
of vortex lines; for stripe solitons, we observed an out-of-phase oscillation of the constituent solitons. Our work paves 
the way for interesting future studies in spin-orbit coupled BECs, e.g., on (line- or ring-) soliton and vortex stability 
and dynamics in higher-dimensional settings. 
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